arXiv: 1505.02480vl [astro-ph.SR] 11 May 2015 


Mon. Not. R. Astron. Soc. 000 ,^-?? (2014) 


Printed 12 May 2015 


(MN IM?eK style file v2.2) 


Survival and Structure of Dusty Vortices in 
Protoplanetary Discs 


Ivo Cmkovic-Rubsamen 1 *, Zhaohuan Zhu 1,2 , and James M. Stone 1 

1 Department of Astrophysical Sciences, 4 Ivy Lone, Peyton Hall, Princeton University, Princeton, NJ 08544 

2 Hubble Fellow 


ABSTRACT 

We have studied the impact of dust feedback on the survival and structure of vor¬ 
tices in protoplanetary discs using 2-D shearing box simulations with Lagrangian dust 
particles. We consider dust with a variety of sizes (stopping time t s = 10" 2 O _1 
10 2 O -1 ), from fully coupled with the gas to the decoupling limit. We find that a vor¬ 
tex is destroyed by dust feedback when the total dust-to-gas mass ratio within the 
vortex is larger than 30-50%, independent of the dust size. The dust distribution can 
still be asymmetric in some cases after the vortex has been destroyed. With smaller 
amounts of dust, a vortex can survive for at least 100 orbits, and the maximum dust 
surface density within the vortex can be more than 100 times larger than the gas 
surface density, potentially facilitating planetesimal formation. On the other hand, in 
these stable vortices, small ( t s < fU 1 ) and large ( t s > f2 _1 ) dust grains concentrate 
differently and affect the gas dynamics in different ways. The distribution of large dust 
is more elongated than that of small dust. Large dust (t s > fl -1 ) concentrates in the 
centre of the vortex and feedback leads to turn-over in vorticity towards the centre, 
forming a quiescent region within an anticyclonic vortex. Such a turn-over is absent if 
the vortex is loaded with small grains. We demonstrate that, in protoplanetary discs 
where both large and small dust grains are present and under the right condition, the 
concentration of large dust towards the vortex centre can lead to a quiescent centre, 
repelling the small dust and forming a small dust ring around the vortex centre. Such 
anticorrelations between small and large dust within vortices may explain the dis¬ 
crepancy between ALMA and near-IR scattered light observations in the asymmetric 
region of transitional discs. 
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1 INTRODUCTION 

One obstacle for planet formation in the core accretion sce¬ 
nario is the growth of solids from millimeter sized grains to 
kilometer planetesimals within the lifetime of protoplane¬ 
tary discs. Various ideas have been proposed, including dust 
collisional coagulation (Weidenschilling 1980, 1995; Blum & 
Wurm 2008 for a review), the gravitational instability for 
dust layers (Safronov & Svjagina 1969; Goldreich & Ward 
1973), the streaming instability (Youdin & Goodman 2005; 
Youdin & Johansen 2007, Johansen & Youdin 2007), and 
particle trapping in large scale disc structures, such as vor¬ 
tices (Barge & Sonmieria 1995; Lyra et al. 2009; Johansen et 
al. 2004; Heng & Kenyon 2010; Meheut, Keppens & Kasse 
2012 ). 

Vortices are often the outcome of the hydrodynamic in- 
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stabilities in rotating flows. In protoplanetary discs, many 
instabilities have been proposed to generate vortices, such 
as the Papaloizou-Pringle Instability (Papaloizou & Pringle 
1984, 1985; Goldreich et al. 1986, Goodman et al. 1987), the 
Rossby wave instability (RWI) (Lovelace et al. 1999; Li et 
al. 2000, 2001; de Val-Borro et al. 2006; Lin & Papaloizou 
2010; Lin 2012a, 2012b, Lovelace & Romanova 2013 for a re¬ 
view), the baroclinic instability (Klahr & Bodenheimer 2003; 
Petersen et al. 2007a, 2007b; Lesur & Papaloizou 2010), Ver¬ 
tical Shear Instability (Urpin & Brandenburg 1998; Nelson 
et al. 2013) instability at buoyancy critical layers (Marcus 
et al. 2013), convection (Godon & Livio 2000), and decaying 
turbulence (Bracco et al. 1999; Shen et al. 2006). 

Two-dimensional simulations on protoplanetary discs 
have shown that anti-cyclonic vortices are long-lived (Godon 
& Livio 1999), and that such vortices can concentrate dust 
particles significantly (Barge & Sommeria 1995; Adams & 
Watkins 1995; Tanga et al. 1996; Godon & Livio 2000), per¬ 
haps thereby facilitating planetesimal formation. 
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However, as dust concentrates in the vortex, it becomes 
more and more critical to consider the dust feedback. When 
the dust-to-gas mass ratio approaches one, it can lead to vor¬ 
tex destruction (Johansen, Andersen & Brandenburg 2004). 
Earlier simulations by Johansen et al. (2004) only evolved 
the vortex for one orbit and the vortex did not reach a steady 
state. In order to concentrate dust significantly, the vortex 
needs to be long-lived. However, when dust is gradually con¬ 
centrated towards the centre of a vortex the feedback be¬ 
comes stronger and the vortex can potentially be destroyed. 
During the preparation of this manuscript, Fu et al. (2014) 
has shown that dust feedback can significantly shorten the 
lifetime of the vortex up to a factor of 10. Raettig et al. 
(2015) have also investigated the stability conditions, but 
that work used convective overstability (Klahr & Hubbard 
2014) to generate vortices. 

In this paper, we carry out a systematic study of dust 
feedback on vortex structure and survival using hydrody- 
namical simulations for up to 100 orbits. The feedback from 
various sized dust with different dust-to-gas ratio is studied. 
In §2, we introduce our method. The results are presented in 
§3. After a short discussion in §4, our conclusions are given 
in §5. 


2 METHOD 

The gas disc was simulated with Athena (Stone et al. 
2008), a higher-order Godunov scheme for hydrodynamics 
and magnetohydrodynamics using the piecewise parabolic 
method (PPM) for spatial reconstruction (Colella & Wood¬ 
ward 1984), and the corner transport upwind (CTU) method 
for multidimensional integration (Klein, Colella & McKee 
1990) (Gardiner & Stone 2005, 2008). Dust particles were 
simulated with the particle integrator in Athena (Bai & 
Stone 2010). Dust particles were implemented as Lagrangian 
particles with the dust-gas coupling drag term, following 


where V; and v 9 denote the velocity vectors for particle i 
and the gas, L is the gravitational force experienced by par¬ 
ticle i, and t s is the stopping time for this particle due to gas 
drag. The dimensionless stopping time T s is defined as 
In this work, we chose the triangular-shaped cloud (TSC) 
interpolation scheme to interpolate the gas properties and 
dust feedback force at the particle position. We use the semi- 
implicit scheme to integrate particle orbits. 

The simulations were conducted using the shearing box 
approximation. The box extends from -0.5 H to 0.5 H with 
256 cells in the x-direction (r-direction) and -H to H with 
512 cells in the y-direction (^-direction), where H is the disc 
scale height (c s /H). We tested this choice of resolution by 
running diagnostic tests at double the resolution (512 by 
1024) and observed no significant change, leading us to be¬ 
lieve that our current resolution was sufficient. We applied 
periodic boundary conditions in the y-direction. To avoid 
wave reflection from x boundaries, densities and velocities 
of ghost zones at x boundaries were fixed at the initial val¬ 
ues (as in Zhu et al. 2014). We set up the Kida vortex by 


initializing the disc structure as 

v x = tt v y/x 

v y = -n v xx , ( 2 ) 

for the region within the ellipse described by Ax and yAx 
as the semi-minor and semi-major axes. In the Kida solu¬ 
tion, Q. v = 3no/(2(x — 1)) (Chavanis 2000). Beyond this 
region, we assumed that the difference between the speed 
given by |2| and the Keplerian speed decreases exponentially 
as exp(-(/ — 1)) where / = (x 2 /Ax 2 + y 2 /( Axx ) 2 ) 1//2 is the 
distance between the vortex centre and the vortex stream¬ 
line in the x direction. In this study we choose Ax = 0.06 
H and x=4- Lyra & Lin (2013) has a good discussion of 
the difference between the Kida solution and the Goodman- 
Narayan-Goldreich solution. We note that for our choice of 
X, the solutions are quite similar. We initially applied a small 
viscosity (v = 10 -5 ) to stabilize the vortex by damping any 
oscillations or sound waves generated by discretizing the in¬ 
compressible Kida solution for use in our compressible grid. 
After 10 orbits when the vortex is stable, the viscosity was 
set to be 0. All runs were conducted to 100 orbits. 

The main set of simulations were conducted to explore 
the parameter phase space given by varying the ratio of solid 
mass to gas mass, e, in the disc and the stopping time of the 
particles. We chose seven e from 0.01 to 1 (e =0.01, 0.0215, 
0.0464, 0.1, 0.215, 0.464, and 1), representative of regions 
of the protoplanetary disc that had already been subject to 
some solid concentration, either by vertical or radial trans¬ 
port. The lower limit e = IIP 2 is the ISM dust-to-gas mass 
ratio. The upper limit e = 1 represents a vortex in a region 
that has already undergone significant solid concentration 
or even streaming instability. 

We also explored the parameter phase space along 
the stopping time dimension. Barge and Sommeria (1995) 
showed that the stopping time is the single parameter that 
characterizes particle concentration in the vortex. The stop¬ 
ping time is dependent on the ratio of mass and the size of 
the particle. Our exploration of stopping time was centered 
on T a = 1. We extended this dimension of the phase space 

IO 4 / 3 

over 3 simulations in each direction from 0.0316 to 
31.6, which captured the limiting behavior in either direc¬ 
tion. 

Thus, we have conducted 7x7 (49) simulations to ex¬ 
plore our phase space, varying the mass ratio and stopping 
time respectively. 

3 RESULTS 

3.1 Initial Stability Assessment 

The long-term stability of vortices for each set of parame¬ 
ters is shown in Fig. [T] The left panel shows the vorticity 
of all runs at 100 orbits. Since all vortices with e > 0.1 are 
destroyed after 100 orbits, we do not show those cases in 
this plot. Throughout the paper, the vorticity is defined as 
uj z — (V x v) z — 1.5f2o, so that the vorticity due to the Ke¬ 
plerian shear has been removed. When the dust to gas mass 
ratio increases, there is a change in the stability of vortices 
near T„ = 1. At high e, only vortices with very large or 
small T s are stable. At very low stopping times the dust 
and gas are totally coupled, while for very large stopping 
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times the dust and gas are uncoupled. Both extremes en¬ 
gender stability. However, particles with T s ss 1 sink to the 
vortex centre fastest (Chavanis 2000) and these vortices are 
the least stable, requiring e ^ 0.0215 for stability. This in¬ 
dicates that for vortices in accretion discs with higher solid 
mass density, stability is heavily constrained. 

Even in cases in which the vortex is destroyed by the 
dust feedback, there is a region of high dust concentration. 
This high dust concentration is produced by the vortex be¬ 
fore it is destroyed by feedback. Thus, even if the vortex is 
a transient phenomenon, it can still lead to very high dust 
concentration. With a small initial dust to gas mass ratio, 
the vortex can survive over 100 orbits independent on the 
size of the particle in the disc. However, the structure of the 
vortex does depend on the size of the particle within the vor¬ 
tex. With small particles ( T a < 1) in the vortex, the vortex 
centre has a minimum vorticity. With big particles (T a >1) 
in the vortex, the vorticity starts to increase towards the 
vortex centre. This difference has important applications on 
dust distribution in vortices, which will be discussed in §4. 

For the smallest particles ( T a = 0.0316) in a stable 
vortex (the uppermost panel of Fig. [I]), we have seen “fin¬ 
gers” at the vortex edge, indicating some instability for well- 
coupled particles. We notice that in a recent paper by Fu et 
al. (2014), a similar phenomenon has been observed. This 
instability may be due to the “heavy core instability” sug¬ 
gested by Chang & Oishi (2010) or parametric instability 
for dust-laden vortices suggested by Railton & Papaloizou 
(2014). 

3.2 Stable Vs. Unstable vortices 

We have studied the time evolution of the vorticity of an 
unstable vortex to understand how decay proceeds. As can 
be seen in Fig. [2] the vortex quickly starts to stretch in the 
y direction until it is strongly deformed. This stretch is a 
consequence of the vortex’s own rotational motion becoming 
slow compared to the shear flow, which is constant through 
the life of the vortex. Once its own vorticity is no longer 
large compared to the shear flow, the vortex is simply pulled 
apart by the shear flow, whereas a stable vortex eventually 
develops a hole as it fills its core with dust but keeps spinning 
relative to the shear flow. For dust with T s > 1, the regions 
of largest spin are at mid-radial distances, as the vortex 
slows down around its dust core. Thus, we conclude that 
even initially stable vortices cannot continue to increase the 
amount of dust they accumulate indefinitely, or they will 
approach the critical value of dust to gas ratio identified in 
§3.3 and destroy themselves. 


3.3 Instability Mechanism Diagnosis 

As particles concentrate within the vortex, the dust feedback 
weakens and eventually destroys the vortex. To quantify the 
weakening of the vortex, Fig.[3]shows the minimum vorticity 
(left panels) and the averaged dust-to-gas mass ratio (right 
panels) within the vortex region. The averaged dust-to-gas 
mass ratio within the vortex is calculated by summing all the 
dust mass within the vortex region and dividing this total 
dust mass by the total gas mass in the same region. The 
vortex region is defined to be an ellipse with a semi-major 


axis of 0.4 H and a semi-minor axis of 0.1 H around the 
vortex centre. In each panel, darker curves represent cases 
with smaller initial dust-to-gas mass ratios. The solid curves 
represent the cases where the vortices are still present at 
100 orbits while the dotted curves represent the cases where 
vortices are destroyed within 100 orbits. Clearly, when the 
vortex is destroyed, the vorticity reaches zero quickly. Again, 
we have observed that the feedback from small dust (e.g., 
T s =0.0316) leads to a smaller vorticity at the vortex centre 
while the feedback from large dust (e.g., T s =31.6) leads to 
a bigger vorticity. 

Fig. [3] also suggests that the survival of the vortex seems 
to depend only on the averaged dust-to-gas mass ratio within 
the vortex, independent of the dust stopping time. For all 
three T a presented in Fig. [3] the vortex survives over 100 or¬ 
bits when the averaged dust-to-gas mass ratio within the 
vortex (right panels) is smaller than 10 -0 ' 5 ~0.3, while 
larger dust-to-gas mass ratios lead to vortex destruction. 

4 DISCUSSION 

4.1 Range of Dust Sizes 

We have identified important differences in the feedback 
caused by small and large dust particles. In a real protoplan¬ 
etary disc, small and large dust coexist in the disc. However, 
the gas will be primarily affected by certain sizes of dust. 
The feedback force from a certain sized dust to the gas is 
the product of the dust abundance and the drag force of 
each dust particle. Generally, a dust species with a higher 
abundance will have a bigger effect on the gas. On the other 
hand, dust with large stopping time will not affect gas sig¬ 
nificantly. Overall, when the gas structure has been altered 
by the feedback from the most important dust species, the 
modified gas structure will affect the distribution of dust at 
other sizes. 

Assuming ISM dust distribution dn(s)/ds ~ s~ 3 ' 5 in 
protoplanetary disks, we know that large dust normally 
dominates most of the dust mass. If this large dust has 
T a > 1, it can make the vorticity increase towards the 
vortex centre (Fig. [T|. With a vorticity close to zero, the 
vortex centre becomes quiescent and small particles can¬ 
not concentrate to the vortex centre. Instead, they will pile 
up around the vortex centre to form a ring structure. Such 
anti-correlation between small and big dust may be able 
to explain the discrepancy of azimuthal structure between 
submm and near-IR scattered light observations for Oph IR 
48 (Follette et al. 2015). 

In order to explore this interesting possibility, we have 
carried out a simulation with particles having a large range 
of sizes from 7),=0.01 to 31.6 with a size distribution of 
dn(s)/ds oc s~ 3 . The total dust-to-gas mass ratio is 0.03. In 
this case, most of the dust mass is in big particles. The gas 
vorticity and dust distribution at 50 orbits are shown in Fig. 
[4] As expected, since large dust concentrates in the vortex 
centre, the dust feedback makes the vorticity approach 0 at 
the vortex centre. This prevents small dust from concentrat¬ 
ing in the vortex centre. Instead, the small particles form a 
ring around the vortex centre. The vortex is gradually weak¬ 
ened by the feedback, and as it becomes more elongated, 
the ring structure of the small particles also elongates. The 
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Figure 1. Left panel: the vorticity of discs with different initial dust-to-gas mass ratios (columns) and different dust stopping times 
(rows) at 100 orbits. The feedback from small dust leads to a vorticity minimum at the vortex centre while the feedback from large 
dust leads to a vorticity turnover at the vortex centre. Right panel: the ratio between dust surface density at 100 orbits and initial dust 
surface density for different cases. The distribution of large dust is more elongated than that of the small dust. 
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Figure 2. Vorticity evolution in a stable (left panel) and unstable (right panel) vortex. The frames are at 5, 10, 20 and 30 orbits, 
descending from top. The stable vortex is from a simulation with mass ratio = 0.0215 and T s = 1, while the unstable vortex has mass 
ratio = 0.0464 and T s = 1. Only the midsection of the shearing box is pictured in each frame for clarity. 


implications of these particle distributions are discussed in 
S4.2. 


4.2 Particle Distribution Within the Vortex 

By comparing ALMA synthetic images based on simula¬ 
tions with ALMA observations towards Oph IRS 48 (van 
der Marel et al. 2014) and HD 142527 (Casassus et al. 2013), 
Zhu & Stone (2014) point out that the vortex centre in these 
observations may be optically thick. Since in regions that are 
optically thick, the intensity is no longer proportional to the 
total dust mass, it is important to know the fraction of the 
vortex region that is optically thick. Thus, we show the dust 
surface density in Fig. [5] The figure indicates that the max¬ 
imum dust density within a stable vortex can reach >100, 
even though this only occurs in a small region around the 
vortex centre. As expected, particles with T a ~1 reach the 
largest concentration at the vortex centre. By assigning an 
opacity to the particles, we know the £d beyond which the 
vortex becomes optically thick. Then with Fig. [5} we can es¬ 
timate how much is hidden in the optically thick region. We 
also note that such high concentration at the vortex centre 


can lead to gravitational instability to form planetesimals or 
even planets (Barge and Sommeria 1995). 


4.3 Vortices with different sizes 

In our main set of simulations, we initialized the vortex with 
a semi-minor axis of 0.06 H. In order to explore how the 
vortex size affects our results, we carried out a similar set 
of simulations with a vortex 5 times bigger (Ax=0.3 H). 
The whole simulation domain was thus adjusted to [-2.5H, 
2.5H]x[-5H, 5H], To save the computational cost, we only 
set the mass ratio (e) to be 0.01, 0.0215, 0.0464 and 0.1, 
and the dust stopping time (T s ) to be 0.0316, 1, 31.6, All 
other parameters, including grid numbers and simulation 
time, are the same as our main set of simulations. We again 
found that vortices with e ^ 0.1 are unstable. Similar to 
Fig. © Fig. © shows that the vortex is destroyed when the 
averaged dust-to-gas mass ratio within the vortex is larger 
than ~0.3. 
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Figure 3. Left panels: the minimum vorticity within the vortex for discs having different sized dust particles. In each panel, darker curves 
represent the cases with smaller initial dust-to-gas mass ratios. The solid curves represent the cases where the vortices are still present 
at 100 orbits while the dotted curves represent the cases where vortices are destroyed within 100 orbits. Right panels: the averaged 
dust-to-gas mass ratio within the vortex region for discs having different sized dust particles. The curves’ color and type have the same 
meaning as the left panel. The vortex is present in cases when the averaged dust-to-gas mass ratio within the vortex is smaller than 0.3. 
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Figure 4. The vorticity (the leftmost panel) and dust distribution (right three panels) from the simulation with particles having a large 
range of sizes from T s = 0.01 to 31.6 with a size distribution of dn(s)/ds oc s ~ 3 . Since large dust concentrates to the vortex centre, the 
dust feedback makes the vorticity increase to 0 at the vortex centre. This makes small particles form a ring around the vortex centre. 
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Figure 6. Similar to Fig. [3] but the vortex is 5 times larger. Similar to small vortices, the vortex is present in cases when the averaged 
dust-to-gas mass ratio within the vortex is smaller than 0.3. 
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Figure 5. The fraction of the vortex region having dust surface 
density larger than E,/ for stable vortices with T s =0.0316. 1, and 
31.6. 


5 CONCLUSIONS 

We have studied dust feedback in vortices in protoplanetary 
discs using 2-D shearing box simulations with Lagrangian 
dust particles. Dust with a variety of sizes (stopping time 
t s = 10~' 2 n _1 - 10 2 fl -1 ), from fully coupled with the gas to 
the decoupling limit, have been considered. 

• We find that the vortex is destroyed by dust feedback 
when the dust-to-gas mass ratio within the vortex is larger 
than 30-50%, independent of the dust size. 


• Even in cases where the vortex is destroyed by dust 
feedback, we find that there is a region of high dust concen¬ 
tration. This concentration is produced by the vortex before 
it is destroyed by feedback. Thus, even if the vortex is a 
transient phenomenon, it can still lead to very high dust 
concentrations. 

• We find that small ( t s < fl -1 ) and large ( t s > D _1 ) 
dust concentrate differently in the vortex and affect the gas 
vortex structure differently. The distribution of large dust 
is more elongated than that of small dust. Large dust (t s > 
QT ) concentrates in the vortex centre and dust feedback 
leads to an turn-over in vorticity towards the vortex centre, 
forming a quiescent centre within the anticyclonic vortex. 
Such a turn-over is absent if the vortex is loaded with small 
dust. 

• We have demonstrated that in protoplanetary discs 
where both small and large dust is present, with the large 
dust comprising most of the dust mass, the concentration 
of large dust towards the vortex centre will lead to a qui¬ 
escent centre, repelling the small dust and forming a small 
dust ring around the vortex centre. Such anticorrelation be¬ 
tween small and large dust within vortices may explain the 
discrepancy between ALMA and near-IR scattered light ob¬ 
servations in the asymmetric region of transitional discs. 

• We find that the dust density can be more than 100 
times larger than the gas density at the vortex centre, po¬ 
tentially facilitating planetesimal and planet formation. 

However, we caution that our 2-D simulations have not 
considered the 3-D structure of the vortex and how verti¬ 
cal transport affects particle distributions. These questions 
of 3-D structure deserve future study. Additionally, we have 
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ignored the radial drift of the particles in our simulations. 
We note that the vortices we study tend to appear at density 
bumps, due to their formation by the Rossby Wave Instabil¬ 
ity for example. These density bumps have much reduced or 
negligible radial drift, so this may limit the magnitude of the 
effect we’re failing to capture. On the other hand, dust radial 
drift may trigger streaming instability which deserves future 
investigation. Finally, the formation of real vortices is driven 
by the presence of a physical instability, which may have sig¬ 
nificant effects on vortex stability and character. Conducting 
simulations which better represent this reality than simply 
introducing vortices as an initial condition could be a fruit¬ 
ful path for future work. Refer to Raettig et al. (2015) for a 
study in this direction. 
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